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Abstract. We investigate bistability and memory effects in a molecular junction weakly coupled to metallic 
leads with the latter being subject to an adiabatic periodic change of the bias voltage. The system is 
described by a simple Anderson-Holstein model and its dynamics is calculated via a master equation 
approach. The controlled electrical switching between the many-body states of the system is achieved 
due to polaron shift and Franck- Condon blockade in the presence of strong electron- vibron interaction. 
Particular emphasis is given to the role played by the excited vibronic states in the bistability and hysteretic 
switching dynamics as a function of the voltage sweeping rates. In general, both the occupation probabilities 
of the vibronic states and the associated vibron energy show hysteretic behaviour for driving frequencies 
in a range set by the minimum and maximum lifetimes of the system. The consequences on the transport 
properties for various driving frequencies and in the limit of DC-bias are also investigated. 

PACS. 85.65.+h - 73.23.-b - 73.40.Gk - 73.63.-b 



1 Introduction 

Quantum switching, bistability and memory effects pro- 
vide potential applications for molecular electronics [TJEl 
[HIS]. Recent scanning-tunneling microscopy (STM) exper- 
iments [5,6,7,8,9] have shown bistability and multista- 
bility of neutral and charged states. Random and con- 
trolled switching of single molecules [TUl[TTj [T2] , as well as 
conformational memory effects [6,9,13,14] have been re- 
cently investigated. Other groups have observed memory 
effects in graphene [T5l[T6l [TT] and carbon nanotubes [T3J 
n [T§1 I5U ] . Motivated by the experimental achievements, sev- 
eral groups f^f^f^fMf2^f26lf27] have attempted to the- 
oretically explain these striking features invoking a strong 
electron- vibron coupling. In Ref. [21 J charge- memory ef- 
fects have been investigated in a polaron-modeled system 
using the equation-of-motion method for the Green's func- 
tions in the strong tunnel coupling regime. Similarly, in 
Ref. [23J these effects are associated with a polaron sys- 
tem treated within a simple mean-field approach. How- 
ever, the hysteresis effects in Ref. [23] may be an arte- 
fact of the mean-field approximation as pointed out by 
Alexandrov and Bratkovsky [58]. In Ref. [53] memory ef- 
fects have been found in a polaron-modeled system tak- 
ing the quantum dot as a d-fold-degenerate energy level 
weakly coupled to the leads and accounting for attractive 
electron-electron interactions. However, here a multiple 
degenerate energy level (d>2) is required. In contrast, in 
Ref. [56] , again the situation of weak coupling to the leads 
but with repulsive electron-electron interaction is consid- 
ered. In this work, bistability, charge-memory effects and 



switching between charged and neutral states of a molecu- 
lar junction have been explained within the framework of 
a polaron model, where an electronic state is coupled to a 
single vibronic mode. These features have been associated 
with the asymmetric voltage drop across the junction and 
the interplay between time scales of voltage sweeping and 
quantum switching rates between metastable states in the 
strong electron- vibron coupling regime. In the weak tun- 
nel coupling limit, a perturbation theory in the tunneling 
amplitude between the molecule and leads is appropri- 
ate to describe electronic transport. In particular, such a 
perturbative treatment is valid if the tunneling-induced 
level width hT is small enough compared to the thermal 
energy k^T. The lowest order in this expansion leads to 
sequential tunneling, which corresponds to the incoherent 
transfer of a single electron from a lead onto the molecule 
or vice versa. Moreover, it is known from transport theory 
that sequential tunneling is dominant as long as the dot 
electrochemical potential (i.e. the difference En — En-i 
between eigenvalues of the many-body Hamiltonian corre- 
sponding to states with particle number differing by unity) 
is located between the Fermi energies of the leads. 

A strong electron-vibron coupling can in turn quali- 
tatively affect the sequential tunneling dynamics [5 ^1301 
[3 H I55 1 l33 l l3l ] ■ For strong coupling, the displacements of 
the potential surfaces for the molecule in a charged or 
neutral configuration are large compared to the quan- 
tum fluctuations of the nuclear configuration in the vibra- 
tional ground state. As a result, the overlap between low- 
lying vibronic states is exponentially small. This leads to 
a low-bias suppression of the sequential transport known 
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as Franck- Condon (FC) blockade, which in turn is respon- 
sible for bistability effects in [26] . 

In this paper we extend and improve the ideas of Ref. [26] 
Specifically, we include the time dependence of the bias 
voltage explicitly, and derive a time-dependent master equa- 
tion for the reduced density matrix of a single level molecule 
coupled to a vibrational mode and weakly coupled to metal- 
lic leads. Moreover, we relax the assumption of fast relax- 
ation of vibrons into their ground states and discuss the 
role played by the vibronic excited states in the switching 
dynamics. As in Ref. [26], we find that controlled electri- 
cal switching between metastable states is achieved due to 
polaron shift and Franck- Condon blockade in the presence 
of strong electron- vibron interaction. Moreover, we find 
that the hysteresis effects can be observed in the switch- 
ing dynamics only if the time scale of variation of the 
external perturbation, T ex , is constrained into a specific 
range set by the minimum, r m i n , and maximum, r max , 
charge lifetimes of the system as a function of the ap- 
plied bias. With A being the dimensionless electron- vibron 
coupling, it holds r m i n ~ r max ~ T -1 e A . Hence, a 
strong electron- vibron coupling (A > 1) is a necessary con- 
dition for the opening of this time scale window and thus 
of hysteresis. Such a large dimensionless electron-vibron 
coupling is not rare in conjugated molecules with soft 
torsional modes (e.g biphenyl with different substituents, 
azobenzene) which have been experimentally proven to 
behave as conformational switches ( [H] , [13] ) . Very large 
reorganization energies (of the order of 1 eV) attributed 
to a polaron effect have also been observed in STM single 
atom switching devices [5]. Also in this case the electron- 
phonon coupling should be large (A > 1) to justify the 
bistability. Outside this range the averaging over multiple 
charging events in the slow driving case or multiple driv- 
ing cycles in the fast case removes the hysteresis. 

The paper is organized as follows: In Section [2] the 
model Hamiltonian of a single level molecule coupled to 
a vibronic mode is introduced. A polaron transformation 
is employed to decouple the electron-vibron interaction 
Hamiltonian and obtain the spectrum of the system. 

In Section [3] we derive equations of motion for the re- 
duced density matrix for the case in which the leads are 
subject to an adiabatic bias sweep. The time-dependent 
master equation is solved in the limit of weak coupling to 
the leads and important time scale relations are derived. 

In Sections EJ [5] and [6j our main results of the memory 
effects are presented and analyzed for a sinusoidal pertur- 
bation of period T ex = 2tt/uj. 

In Section [4] the lifetimes of the many-body states of the 
system are calculated. We show that, for the case of asym- 
metric voltage drop across the junction, at small bias volt- 
ages a bistable configuration is achieved which plays a 
significant role in the hysteretic dynamics of the system. 
Bistability can involve also vibronic excited states of the 
system. 

In Sections [5] and [6] we give an explanation of the hys- 
teretic behavior of the system in terms of characteristic 
time scales, in particular, the interplay between the time 
scale T ex of variation of the external perturbation and of 



the dynamics of the system set by r sw i tc h ~ Tmin ~ 
In Section [5] focus is on the regime uj ~ T while in Sec- 
tion [6] is uj <C T. In the latter case the features observed 
in Ref. [56] can be successfully reproduced. 

In Section[7J the consequences on the transport proper- 
ties in the DC-limit are presented as a special case. Finally, 
we conclude in Section [8] 



2 Model Hamiltonian 

We consider a simple Anderson-Holstein model where the 
Hamiltonian of the central system is described as 

H sys = H mo \ + H v + #e-v, (1) 

where H mo \ represents a spinless single molecular level 
modeled by the Hamiltonian 

H mo i = (e + eV e )$d, (2) 

where w{d) is the creation (annihilation) operator of an 
electron on the molecule and eo is the energy of the molec- 
ular level, and V g accounts for an externally applied gate 
voltage. For simplicity we assume a spinless state describ- 
ing the molecular level with strong Coulomb interaction 
where only one excess electron is taken into account. The 
spin degeneracy would not qualitatively change the results 
of the paper. The vibron Hamiltonian can be written as 

H v = ficjQ (^a + ^ , (3) 

where a^(a) creates (annihilates) a vibron with energy 
hjjQ. Finally, the electron-vibron interaction Hamiltonian 
is expressed as 

ff^ v = 0<?d(at + a), (4) 

where g is a coupling constant. 
2.1 Polaron transformation 

In order to decouple the electron-vibron interaction Hamil- 
tonian, we apply the canonical polaron unitary transfor- 
mation [35]. Explicitly, we set H S ys — e S ' H sys e~ s , where 

S = \$d(tf - a), (5) 

with A = as the dimensionless coupling constant. The 
transformed form of the electron operator is 

d = dX, (6) 

where X = exp [—A (aJ — a)] . In a similar way, the vibron 
operator is transformed as 

a = a - \d)d. (7) 
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Now the transformed form of the system Hamiltonian reads 



^sys = eSd + huuo ( o) 



1 

■a + - 



(8) 



where e = eo-\-eV g — -fj^ is the polaron energy with polaron 

2 

shift Er, = S~r. The polaron eigenstates of the system are 



hojQ 



\n,m)i := 



n,m), 



(9) 



where n denotes the number of electrons on the molecular 
quantum dot, while the quantum number m characterizes 
a vibrational excitation induced by the electron transfer 
to or from the dot. 



3 Sequential tunneling 

We analyze the transport properties of the system in the 
limit of weak coupling to the leads. The Hamiltonian of 
the full system is expressed as 



H(t) = H sys + Ht 



(10) 



where a = s,d, denotes the source and the drain contacts, 
respectively. The tunneling Hamiltonian is given by 

H T = J2 ta (cLd + d f Ca«) , (11) 



where c^ aK (c aK ) creates (annihilates) an electron in lead a. 
The coupling between molecule and leads is parametrized 
by the tunneling matrix elements t s and t d . Here, we con- 
sider the weak coupling regime so that the energy broad- 
ening hT of molecular levels due to Ht is small, i.e., hT <C 
hjOQ^k^T^ and a perturbative treatment for Ht in the 
framework of rate equations is appropriate. For simplicity, 
we assume that the tunneling amplitude t s / d of lead s/d 
is real and independent of the momentum Tin of the lead 
state. In addition, we consider a symmetric device with 
t s = td' Finally, the time dependent lead Hamiltonian is 
described by 



H a (t) = ^2[e K + AfjL a {t)]i 



(12) 



The above equation describes the lead Hamiltonian of non- 
interacting electrons with dispersion relation e K . The time- 
varying chemical potential Afi a (t) of lead a depends on 
the applied bias voltage, and yields a ^-independent shift 
of all the single-particle levels. 



3.1 Time dependent master equations for the reduced 
density matrix 

In this section, we briefly derive the equation of motion 
for the reduced density matrix (RDM) of the molecular 



junction accounting for the time-dependence, Eq. ([T2|h 
of the lead Hamiltonian H a (t). We restrict to the lowest 
nonvanishing order in the tunneling Hamiltonian. Never- 
theless, due to the explicit time dependance in the leads 
Hamiltonian, this work represents an extension of previ- 
ous studies on similar systems (see e.g., Refs. [33 J3 ^[3fH 
[STIIMIMIira The method is 

based on the well known Liouville equation for the time 
evolution of the density matrix of the full system consist- 
ing of the leads and the generic quantum dot. To describe 
the electronic transport through the molecule, we solve 
the Liouville equation 



in* t — -Lrieads 

at 



H^J^t) 



(13) 



for the reduced density matrix p re d(t) = Tri ea d s {p(t)} in 
the interaction picture, where the trace over the leads de- 
grees of freedom is taken. In the above equation, H^(t) is 
the tunneling Hamiltonian in the interaction picture to be 
calculated as below: 



H^it) = J2t a \clJ(t)e^ t+ ^ +h.c. 



(14) 



where ( a (t) — f* Ap a {t')dt' . We make the following ap- 
proximations to solve the above equation: (i) The leads 
are considered as reservoirs of noninteracting electrons in 
adiabatic thermal equilibrium. Note that this implies that 
the time scale of variation of the external perturbation 
has to be large compared to the relaxation time scale of 
the reservoirs (cf. Eq. (fl9|) below). We assume the cou- 
pling between system and reservoirs has been switched on 
at time t = to and consider a factorized initial condition. 
Thus at times t > to it holds p 1 (t) = pl ys (t) (8) p s pd + 

0(t - t )O(H T ) := pi ys (t) /wis + 0(t - t )O(H T ), where 
the correction in the tunnelling Hamiltonian drops in the 
second order master equation (see Eq. ([T6]) ). Here p s / d = 
_i_ e -!3(H s/d (t)-fi s/d (t)N s/d ) d eno t es the thermal equilibrium 



grandcanonical distribution of lead s/d, Z s / d is the par- 
tition function, f3 the inverse of the thermal energy, 7V s / d 
the electron number operator, and p s /d(t) — Mo + A/^s/dW 
is the time dependent chemical potential of lead s/d which 
depends on the applied bias voltage. Note that the levels 
shift is taken into account by the time-dependent pertur- 
bation A/i s / d (£), while the change in chemical potential is 
taken into account accordingly via the chemical potential 
Ps/d(t) so that the net positive or negative charge accu- 
mulation in the leads is avoided. Conventionally, we take 
the molecular energy levels as a fixed reference and let 
the bias voltage drop across the source and drain contacts 
through the Fermi energies as [52] 



Ps(t) = /JL + (1 - r])eV h (t), 
fi d {t) = po ~ rjeVb(t), 



(15) 



where < rj < 1 describes the symmetry of the volt- 
age drop across the junction. Specifically, r] = corre- 
sponds to the most asymmetric situation, while rj = 1/2 
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represents the symmetric case. In addition, we consider a 
sinusoidally- varying bias voltage, i.e., Vb(t) = Vbsin(o;t), 
where uo is the frequency of the driving field, (ii) Since 
we assume weak coupling of the molecule to the leads, we 
treat the effects of Ht perturbatively up to second order. 
Accounting for the time-evolution as in Eq. ([T4|) of the 
leads creation/annihilation operators, we find: 



x e *Mt-0+C«tt-C(0] + [i-f a {£k - /*)] 
xdt W d(t')^ ed ( t ') e -i^(t-t')+a(t)-C»(t')] 

-[l-f a (e K -^)]d(t)pi ed (t')d\t') 

x ^Wd^e-te^-*^®-*^ +h.c|. 

(16) 

In the derivation of the above equation we have used the 
relation: Tri eads {c^c^^psPd} = <W<W/ (s K - /io), 
where / (e K — po) is the Fermi function, and the cyclic 
property of the trace. By summing over n we obtain the 
generalized master equation (GME) for the reduced den- 
sity matrix in the form 

xp4 d (Oe*^W-^(*')]+F a (i-t',- Mo ) 

-^(t-t'.-Aio)^*)^^^) 
xe*K-W-^(*')]-F*(t-t',/io)^(t) 

x /3ied(*')d(*')e~* KB( * ) ~ < " ( *' )1 + h.cj, (17) 

where the correlation function F a (t — t', fio) of lead a [see 
Appendix |X] has, in the wide band limit, the following 
form: 



F a (t-t',p ) 



■■nhD a e i! 2(*-t') 



Hj3 sinh \n 



(t-t') i [' 
hp J J 



(18) 



which decays with the time difference t — t' approximately 



. (*-*') 

h(3 



on the time scale — . Here is the 



as exp 

density of states of lead a at the Fermi level, (iii) Since 
we are interested in the long-term dynamical behavior of 
the system, we set to —> — oo in Eq. ([T7|) . Furthermore, we 
replace t' by t — t" . We then apply the Markov approxima- 
tion, where the time evolution of p^ ed is taken only local 
in time, meaning we approximate p re d(t — t") ~ p r ed(t) in 
Eq. ([T7|). In general the condition of time locality requires 
that [50J 



Here we defined from Eq. ([TT]) together with Eq. ([T8)h 
T a = jr-\t a \ 2 D a as the bare transfer rates and HT = 
hT a as the tunneling-induced level width. Notice that 
the validity of the Markov approximation, justified in this 
case, is crucially depending by the order of the current cu- 
mulant and the order of the perturbation expansion in the 
tunnelling coupling [5T] , Finally, the condition of adiabatic 
driving Eq. ([T9|) allows to approximate ( a (t) —( a (t — t") = 
Ap a (t)t" . Taking into account these simplifications, the 
generalized master equation (GME) for the reduced den- 
sity matrix acquires the form 

tied® = - X) ^ /°° - 

x pl ed (t) + F[t", -p a (t)}d\t)d(t - t")pl ed (t) 

- F*[t", -p a (t)]d(t)pi ed (t)d\t - t") 

- F*[t\ fi a (t)]d\t)pi ed (t)d(t - t") + h.c.l, (20) 



(19) 



where F[t" , p a (t)} = F a (t" , po)e^ a{ ~ t)t " . Since the eigen- 
states |n, m) 1 of H sy s are known, it is convenient to calcu- 
late the time evolution of p^ ed in this basis. For a generic 
quantum dot system, this projection yields a set of differ- 
ential equations coupling diagonal (populations) and off- 
diagonal (coherences) components of the RDM. For the 
simple Anderson-Holstein model Eq. (pp) coherences and 
populations are, however, decoupled. In the sequential- 
tunneling regime, the master equation for the occupation 
probabilities P™ = i(n, m\p re d\n, m)i of finding the sys- 
tem in one of the polaron eigenstates assumes the form 

Pn = E c'-Twc' - E r "'wc- (2i) 



where the inequality F <C ooo ensures the applicability of 
the secular approximation, i.e., the separation between the 
dynamics of populations and coherences. In the numeri- 
cal treatment of these equations we truncate the phonon 
space. Convergence is reached already with 40 excitations. 
In Eq. (|2Tj) the coefficient r™,^ 71 denotes the transition 
rate from \n' ', m')i into the many body state |n, m)i, while 
r^!^ 1 describes the transition rate out of the state |n, m)\ 
to |n', m')i. Taking into account all possible single-electron- 
tunneling processes, we obtain the incoming and outgoing 
tunneling rates, in the wide band limit, as 

rS£ m '(f) = E r * F ™m'f + [e + ^0 K - m) - fl a (t)] 

a 

= E r Ki'W' (22) 



Tf^ m (t) = E r a F mm ,f~ [e + Hw (m' - m) - /*„(*)] 

a. 

= E r ™'i-oW, (23) 
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where the terms describing sequential tunneling from and 
to the lead a are proportional to the Fermi functions 
f+(x - iicx) = f(x - l^cx) and f~(x - fi a ) = 1 - f(x - 
fi a ), respectively. Notice that the integrations over en- 
ergy and time introduce the explicit time dependance in 
the Fermi functions. The factor F mm / = \{m\X\m')\ 2 is 
the Franck- Condon matrix element which can be calcu- 
lated, with X defined in Section [27L] explicitly using Ap- 
pendix The sum rules ]T m F mm / £) m , F mm / 1 
are well satisfied because of the completeness of each vi- 
brational basis set |0,m) and \l,m f )i. This factor de- 
scribes the wave-function overlap between the vibronic 
states participating in the particular transition. It contains 
essential information about the quantum mechanics of the 
molecule and significantly influences the transport prop- 
erties of the single- molecule junction. Within the rate- 
equation approach, the (particle) current through lead a 
is determined by 

i a {t) = J2 (c/Ti'W W) - r£iZS(t)ir'(tj) (24) 

mm' 

and it is in general time dependent. Moreover, differently 
from the stationary case, in general ^ — Iii(t). The 

charge is though not accumulating on the dot since, for 
the average quantities 



^,av lim 



pt+T ey 
Jt 



dt f I a (t f ) 



(25) 



it holds iL,av = —Ir,&v, as it can be easily proved consid- 
ering that the average charge on the dot oscillates with 
the same period T ex of the driving bias. Finally, in the 
DC limit uj —> the relation = holds as the 

fully adiabatic driving allows to reach the quasi-stationary 
limit at all times. 



4 Lifetimes and bistability of states 

In this section, we show that when the bias voltage drop is 
asymmetric across the junction, upon sweeping the bias, 
one can tune the lifetime of the neutral and charged states 
to achieve a bistable system. The lifetime of a state is 
obtained by calculating the switching rate of that state. 
The lifetime r nm of a generic quantum state |n,ra)i is 
given by the sum of the rates of all possible processes 
which depopulate this state, i.e., 



(26) 



and it defines, at least on a relative scale, the stability of 
the state |n, m)\. Thus, at finite bias voltage, the inverse 
lifetime of the O-particle mth vibronic state is given by the 
relation 



'Ora 



r «^mm'/ + [e + huo (m - m) - fj, a ] . (27) 



In a similar way, the inverse lifetime of the 1-particle and 
mth vibronic state is expressed as 

= T ocFmm'f~ [e + hujQ (m - rri) - /i a ] . (28) 

a,m f 



A consequence of Eqs. (|27j) and ([28]) is that, due to the 
characteristic features of the Franck- Condon matrix ele- 
ments, in the strong electron- vibron coupling regime, the 
tunneling with small changes in m — m' is suppressed ex- 
ponentially. Hence only some selected vibronic states con- 
tribute to the tunneling process. However, tunneling also 
depends on the bias voltage and temperature through the 
Fermi function. To proceed further, let us focus first on 
the lifetime of the 0- and 1-particle ground states for the 
case of fully asymmetric coupling of the bias voltage to 
the leads, i.e., r] = 0: 



W = Yl n — { r ^ + ( £ + m '^o - Mo - eV h ) 

771 . 

m 

+ r d /+ (s + m'tuvo- no)}, (29) 



e -X 2 ^2m f 

Tiq 1 = ^2 — — { r s/ _ (e - rn'hujQ - /i - eV h ) 



+ r d / (e - m'hwQ - /io)}- 



(30) 



One can see from Eq. ([29]) that if in the considered param- 
eters range is e + m'hujo ^> Mo, i-e., / (e + m'huo — Mo) — ^ 
0, then the second term in the bracket is negligible. The 
first term is nonzero at large positive bias, while at large 
negative bias it remains negligible. In a similar way one 
can analyze the behavior of r^ 1 in which the first term on 
the r.h.s. of Eq. ([30]) will be dominating at large negative 
bias. In order to understand the mechanism of this pro- 
cess the energy-level scheme for the relevant transitions in 
a coordinate system given by the particle number TV and 
the grandcanonical energy E — hqN shown in Figure [U We 
choose V g = and /io = 0. Moreover, the polaron energy 
levels are at resonance with the 0-particle states for our 
chosen set of parameters: we set e p = eo and hence e = 0. 
Then the only transitions allowed at zero bias are ground 
state o ground state transitions. At finite bias also tran- 
sitions involving excited vibronic states become allowed. 
In particular, at Vh = it follows from Eqs. ([29]) , and ([30]) 
that 



TooH^b = o) = rfoHVb = o) = e - A2 (r s + r d )/2, (31) 

while at |Vb| — > oo it holds 

riioiVb -> °°) = T w(Vi> -> -°°) 

= r s + ^ 2 ^r s ^r-i, (32) 



whereas 



Too 1 ^ -»■ -00) = rfo^Vi, y 00) = ^e" A2 = r," 1 



a,m f 



(33) 
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Fig. 1. (Color online), (a) Energy- level scheme for the rele- 
vant transitions in a coordinate system given by the particle 
number TV and grandcanonical energy E — /noN at Vh = 0. 
The red lines represent the transitions threshold, where the 
thickness of each transition line gives the strength of the tran- 
sition. The polaron energy levels are aligned with the 0-particle 
states for our chosen set of parameters (fio = 0, V g — 0, So — 
25/io;o, A = 5) yielding the polaron shift e p = Sq. (b) Inverse 
lifetimes (r n or) _1 on logarithmic scale as a function of nor- 
malized bias voltage eVh/fowo. The red thick line represents 
the inverse lifetime of the 1-particle ground state, while the 
thin blue line refers to the 0-particle ground state. 




eV b /fiu) 
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Fig. 2. (Color online). Inverse lifetime (r nm r) _1 as a function 
of normalized bias voltage eVb/fowo for (a) vibronic ground 
states, (b) first excited states, (c) second excited states, (d) 
third excited states, (e) fourth excited states, (f) fifth excited 
states when V g = 0. The blue thin line represents the inverse 
lifetime of the 0-particle state (n = 0), while the thick dashed 
red line refers to the 1-particle state (n — 1). The asymmetry 
parameter is r\ = and we fix the zero of the energy at the 
leads chemical potential at zero bias: fio — 0. The energy of 
the molecular level is eo = 25huuo . The electron- vibron coupling 
constant is A = 5 yielding a polaron shift s p = sq. Finally, the 
thermal energy is k&T = 0.2/kJo, the frequency of the driving 
field is u = 0.002cj , and T s = T d = 0.006cj . 
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Fig. 3. (Color online). Inverse lifetime (r nm r) _1 as a function 
of normalized bias voltage eVh/hxiQ for (a) vibronic ground 
states, (b) first excited states, (c) fifth excited states when 
eV g /huJo = 8, while (d) shows vibronic ground states, (e) first 
excited states, and (f) fifth excited states, when eV g /huuo = —8. 
The remaining parameters are the same as used in Figure [2] 



In practice the asymptotic behaviors are already reached 
at e\ Vh \/fiwo ~ 2 A 2 as observed in Figure HJb). Note that 
T max and T m i n set the maximum and minimum achiev- 
able lifetimes which, due to r max /r m i n ~ e A , can differ 
by several orders of magnitude for A > 1. Note also that 
near zero bias the lifetimes are so long that the system 
never likes to charge or discharge and a bistable situation 
is reached. A selective switching, however, can occur upon 
sweeping the bias voltage. Hence r m i n also sets the time 
scale for switching: r min - T switc h. 

Analogously, we can explain the behavior of the life- 
times of the excited states (see Figure [2j). It follows that in 
the considered parameters range, in general, the 0-particle 



vibronic states are stable at large enough negative bias 
voltage, while the 1-particle vibronic states are stable at 
large positive bias. There is, however, an interval of bias 
voltage, the so-called bistable region, where both states 
|l,m')i and |0,ra)i are metastable for not too large m 
and ra/ , as shown in Figure [2j Moreover, m steps are ob- 
served in the inverse lifetimes (see Figures [2fb-f)) 
because for certain values of the coupling constant A some 
of the FC factors F mrn r vanish or are exponentially small 
such that the additional channels opened upon increasing 
the bias voltage do not have pronounced contribution. For 
instance, the FC factor for the first excited vibronic state 
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can be described as 



5 Quantum switching and hysteresis 



_ \2(m'-l) 

F w =e- A2 ^^(m'-A 2 ) 2 , 



m'\ 



(34) 



which vanishes for m' = A 2 . That is why a plateau around 
eVb/hoJo = 25 in Figure [2Kb) is observed for our chosen 
parameters. Analogously, using Eq. ([5T]h one can find (cf. 
Appendix [D]) that i^m 7 has two minima at 



1 + 2A 2 + vTTlA 2 1 + 2A 2 - vTTlA 2 
mi = , ra 2 = . 

(35) 

Hence two plateau can be observed (see Figure[2Kc)) around 
eVh/huJo = 20 and eVh/huo = 31. Similar arguments can 
be extended to explain the steps in the inverse lifetimes of 
higher excited states. This also implies that the bias win- 
dow for bistability shrinks for excited states and even dis- 
appears for large enough m. It follows that the major con- 
tribution in bistability is coming from low excited vibronic 
states. Note that the bistability of the many body states is 
crucial for the hysteresis and hence memory effects which 
is discussed in the next section. Finally, a closer inspection 
of Figure [2] reveals that the minimum of the inverse life- 
time increases with the vibronic quantum number m. This 
effect can be understood easily by analyzing the minimum 
of the inverse lifetime for each particle state. For exam- 
ple the minimum of the inverse lifetime for the 0-particle 
vibronic ground state is, cf Eq. (|33]h whereas for the 0- 
particle vibronic first excited state is 



'01 



1 (F b ^-oo) = ^(l + A 4 )e- 



(36) 



From Eqs. ([33]) and ([36]h one can conclude that T 00 1 (Vb — >• 
— oo ) < TQi(Vb — » — oo ). A similar explanation can be ex- 
tended to the higher excited states. For gate voltages such 
that eV g > 0, the 1-particle vibronic excited states are be- 
coming unstable faster than the O-particle states (see Fig- 
ure EKa)-(c)), while for large negative gate (eV g < 0), the 
O-particle states are getting unstable fast (see Figure [3^d)- 
(f)). In order to explain this effect, we analyze the shift of 
the inverse lifetime of the O-particle vibronic first excited 
state, Tq^ 1 , as follows: 

The maximum of the inverse lifetime for V g =^ is 
Tbi^Vb -> oo) = T s + r d ^F lm /(eVg + fio;o(m - 1)), 

m 

(37) 

whereas the minimum is given by 
roi(V h ->■ -to) = T d J2 Fimf(eV g + hw (m - 1)). (38) 

m 

Eqs. ([37|) and ([38|) imply that both minimum and maxi- 
mum of Tq^ 1 shift by an equal amount and the condition 
of the bistability region can be tuned by setting V g . 



Neutral and charged (polaron) states correspond to dif- 
ferent potential energy surfaces and transitions between 
low-lying vibronic states are strongly suppressed in the 
presence of strong electron-vibron interaction. This leads 
to the bistability of the system. Upon applying an exter- 
nal voltage, one can change the state of this bistable sys- 
tem obtaining under specific conditions hysteretic charge- 
voltage and current-voltage curves. Here it is crucial to 
point out that only if the time scale of variation of the ex- 
ternal perturbation is shorter than the maximum lifetime 
but longer than the minimum lifetime of the system hys- 
teresis can be observed, i.e., r min ~ T switc h < T ex < r max . 
Due to r max > T ex , the system stays in the stable state 
during the sweeping until the sign of the perturbation 
changes, the former stable state becomes unstable and, 
due to T ex < T m i n , a switching to the new stable state can 
occur. In this section we now consider the situation when 
uo ~ r, i.e., T ex ~ r sw itch while in Section [6] the regime 
uj <C T, i.e., T ex > r sw itch is addressed. 




100 150 200 250 

t[i/r] 
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Fig. 4. (Color online), (a)-(b) Occupation probabilities Po and 
Pi of the 0- and 1-particle electronic states as a function of nor- 
malized time dependent bias voltage eVb/huJo, (c) population 
of the O-particle configuration as a function of time; (d) nor- 
malized bias voltage as a function of time. The parameters are 
the same as used in Figure [2] 



In Figures H] and [5] we present the populations of the 
electronic states, P n = P^ 1 , as well as of the vibronic 
states, P m = ^2 n P™, respectively. Specifically, in Fig- 
ure HKa)-(b), we have plotted the populations of the 0- 
and 1-particle electronic states as a function of normal- 
ized bias voltage, where hysteresis loops can be seen. In 
Figure HKc), instead, we have shown the population of the 
O-particle electronic state as a function of time. The lat- 
ter can be used to determine the time T sw itch of switching 
between the neutral and charged states. In a similar way, 



8 



Andrea Donarini, Abdullah Yar, and Milena Grifoni: Vibration induced memory effects and switching in driven. 



the sweeping time T ex of the bias voltage can be calcu- 
lated using Figure Hfd). By comparison of these two time 
scales, it is apparent that the switching time is of the 
same order as the sweeping time and much shorter than 
the lifetime in the bistable region (see Figure [T]). The rela- 
tion T sw itch ~ ^ex also explains why the switching between 
the neutral and charged state is on average never complete 
(Po oscillates between 0.2 and 0.8). 
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Fig. 6. (Color online). Plots of the population P™ as a func- 
tion of normalized time dependent voltage eVb/fowo for the 
0-particle vibronic (a) ground state, (b) first excited state, (c) 
fifth excited state, and for the 1-particle vibronic (d) ground 
state, (e) first excited state, (f) fifth excited state. The param- 
eters are the same as used in Figure [2] 
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Fig. 5. (Color online). Populations P m of the vibronic states as 
a function of normalized time dependent bias voltage eVh / fowo 
for (a) ground state, (b) first excited state, (c) second excited 
state, and (d) fifth excited state. The parameters are the same 
as used in Figured 



5.1 I — V characteristics 

The hysteretic behavior of the bistable system is also re- 
flected in the current as a function of normalized bias (see 
Figure [7]) where a hysteresis loop (single loop) is observed 
in the current calculated both at the left and the right 
lead. Interestingly, the left and the right currents differ by 
more than a sign, in contrast to the stationary case. This 
behavior is understandable again in terms of relaxation 
time scales. In fact, for voltages |V&| outside the bistable 
region the system relaxes to the stationary regime on a 
time scale r sw i tc h- Though, since the driving time T ex has 
the same order of magnitude, the stationary regime can- 
not be reached. Yet, no net charge accumulation occurs 
since 7 L>av = -flav- 



in Figure [5] the populations of the vibronic states as a 
function of the normalized bias voltage are shown, while 
in Figure [6] the populations of the different vibronic states 
resolved for different charges have been plotted. Clearly 
not only the vibronic ground states (which were consid- 
ered in Ref. [26]) show hysteretic behavior but the vibronic 
excited states also exhibit these interesting features. Fur- 
thermore, inspection of these figures reveals that even af- 
ter relaxation on the stable limit cycle, the vibronic ex- 
cited states are highly populated in the non-stationary 
case in contrast to the stationary case uo — >• (see e.g., 
Figures [15] and [T7|) where the population of the excited 
states is strongly suppressed. Finally, while the general 
trend is a reduction of the population, the higher the ex- 
citation and the populations are negligible for m « 40, an 
interesting behaviour can be recognized in the form of the 
limit cycles. Namely, upon sweeping the bias we find that, 
for m ^> 8 the probability grows at large biases, it stays 
essentially constant for m « 8 and it decreases at larger 
biases for m < 8. The interpretation of this behaviour is 
still unclear to us. All these observation confirm, though, 
that it is natural to take into account the vibronic excited 
states in the dynamics of the system. 




-60 -30 30 60 -60 -30 30 60 



eVb/hu) eVb/hcJo 

Fig. 7. (Color online). Time dependent current as a function 
of normalized voltage for (a) left lead, (b) right lead. The pa- 
rameters are the same as used in Figure [2] 

In Figure [8] we plot the left time dependent current 
as a function of the normalized bias for different values 
of the electron- vibron coupling constant. An inspection of 
this figure reveals that the width of the hysteresis loop de- 
creases and shifts from zero bias upon decreasing the cou- 
pling constant A. This feature can be understood by ob- 
serving that for A ^ 5 the polaron shift e p does not longer 
compensate the energy of the molecular level £o? and hence 
the polaron energy e ^ 0. In other words, the system is 
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no longer behaving symmetrically upon exchange of the 
sign of the bias voltage. If we consider e.g. the case A = 1 
is, for V g = no — 0, s/hwo = 24. In turn this implies that 
^(^b = 0) - and rfo^Vb = 0) ~ T s + r d , i.e., the re- 
gion around zero bias is no longer bistable as for the case 
A = 5. Hence the dot is preferably empty at zero bias. 
Switching however can be reached upon increasing Vh in 
the region around eVh ~ e. Overall however the bistabil- 
ity region has shrunk. Similar considerations apply to the 
other considered values of A. 



0.2 

l_ 

0^ o 
_ l 

-0.1 



-0.2 



(a) 






A=4 



-60 



-30 30 

eVb/ticoo 



60 




-30 30 

eVb/ficoo 




-30 30 

eVb/fiooo 



30 30 

eVb/fiooo 



Fig. 8. (Color online). Time dependent current for the left lead 
as a function of the normalized voltage for coupling constants 
(a) A = 4, (b) A = 3, (c) A = 2, and (d) A = 1. The remaining 
parameters are the same as used in Figure [2] 



The vibron energy associated with the 0-particle state 
is determined by the relation 



^v,o — Tr sys 



(40) 



with po = p re d|0,ra)ii(0, m\. In Figure [51(b), the normal- 
ized vibronic energy as a function of normalized bias volt- 
age for the 0-particle configuration has been plotted. The 
hysteresis loop resembles that of Figure H]^a) implying a 
direct correlation between the vibronic energy and the 
population of the neutral state i.e., the more the neu- 
tral state is occupied the higher is the associated vibronic 
energy. Qualitatively the result can be explained as fol- 
lows: transitions from the charged to the neutral states 
are predominantly involving low energy charged states and 
highly excited neutral states. Due to energy conservation 
and asymmetric bias drop these transitions are confined to 
the large negative biases where the highly excited neutral 
states show also a long life time. This situation remains 
roughly unchanged during the up sweep of the bias until 
the symmetric condition is obtained at high positive bias 
and the charged excited states are maximally populated. 
Finally, the bistability around zero bias explains the hys- 
teresis. The analytical expression for the vibronic energy 




eV b /fioo 



eV b /Ficj 
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Fig. 9. (Color online), (a) Total vibron energy as a function 
of the time dependent bias voltage, (b) Vibron energy for the 
0-particle, and (c) for the 1-particle configuration only. Param- 
eters are the same as used in Figure [2] 



5.2 Vibron energy 

In this section, we illustrate the role played by the vibronic 
energy in the hysteretic behavior of the system. The vi- 
bron energy of the whole system can be expressed as 




where the trace is taken over the system degrees of free- 
dom. The normalized vibronic energy as a function of nor- 
malized bias voltage is depicted in Figure [9fa), where hys- 
teretic loops are also observed. The value of the vibronic 
energy, together with the observation that the probabil- 
ity distribution is relatively flat over the excitations (see 
Fig. [6j) ensures that, depending on the bias, between 10 
and 20 vibronic excited states are considerably populated. 
Further insight in the dynamics of the system is obtained 
by considering the correlation between the vibronic energy 
and the charge occupation. 



of the 1-particle state is given by 




with pi = p re d|l, m)n(l, m|. The normalized average vi- 
bron energy as a function of normalized bias voltage for 
the 1-particle configuration is sketched in Figure [9^c), 
where we can observe a hysteresis loop resembling that 
of Figure [3Jb). 

In conclusion, the vibron energies also show hysteretic 
behavior, in analogy to the population- voltage and current- 
voltage curves, in the non- stationary limit. 

6 Testing lower driving frequencies 

When lowering the driving frequency uo (uj <C T) of the ex- 
ternal perturbation, we choose w = 2x 10 _6 o;o, our model 
displays features similar to those presented in Ref. |26| . 
In more detail, we show the population of the electronic 
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states as a function of normalized bias and time in Fig- 
ure [10^ a) -(b), Figure ITPlf c), respectively, whereas in Fig- 
ure [TUT d) the normalized bias as a function of time is 
shown. In this case the population- voltage curve is slightly 
different from Figure H] because the transition between 
and 1 occurs more abruptly as a function of and it is 
complete. Indeed, for the parameter chosen in Figure [TOl is 
e = and r~;[ x ~ cj min <Cw<r^ r S witch- ^ n °ther words 
the frequency is small compared to the charge/discharge 
rate. The system thus follows adiabatically the changes of 
the bias voltage and only switches at those values of the 
bias where r n o ~ T sw i tc h- The time-dependent left current 
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Fig. 10. (Color online), (a)-(b) Population of the 0- and 1- 
particle electronic states as a function of the bias voltage, (c) 
population of the 0-particle electronic state as a function of 
time, and (d) normalized bias voltage as a function of time. 



The frequency of the driving field is uj — 2 x 10 
other parameters are the same as used in Figure [2] 



The 



as a function of normalized bias is shown in Figure [TTt a) 
giving two loops, one for positive bias sweeping and the 
other for negative sweeping. The right current is shown 
in Figure [TTTb). Due to the extremely low frequency the 
currents substantially fulfill the quasi-stationary relation 
lL,(f) — —Iii(t) associated to a fully adiabatic regime. 
In Figure HH we present the populations of the vibronic 
states and hysteretic loops are visible. Vibronic states with 
quantum numbers up to A all display nonvanishing popu- 
lations, much less than in the case T ex w T sw i tc h- 



7 The DC-case (uj -> 0) 

In this section, we consider the limit (uj — >• 0) of DC- 
bias as a special case of the master equation presented in 
the previous section and compare the results. Even if the 
system still exhibits the bistable properties discussed in 
Section H] (they are in fact not related to the sweeping time 
of the bias) the hysteretic behavior cannot be observed 
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Fig. 11. (Color online), (a) Plot of the time dependent current 
as a function of normalized voltage, (b) current averaged over 
one driving period as a function of normalized voltage. The 
value of gate voltage is V g = and the frequency of the driving 



field is uj = 2 x 10 
used in Figure [2] 



^(jJq. The other parameters are the same as 
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Fig. 12. (Color online). Plots of populations of vibronic (a) 
ground state, (b) first excited state, (c) second excited state, 
and (d) fifth excited state. The frequency of the driving field 
is uj = 2 x 10 _6 cjo- The other parameters are the same as used 
in Figured] 



anymore. In Figure [13l we present the population of the 
electronic states for gate voltage V g = 0. At large negative 
bias the system is empty, while at large positive bias it is 
charged. The system makes transitions from the 0- to 1- 
particle state near zero bias. Analogously, in Figure [I4j the 
population of electronic states as a function of normalized 
bias is depicted for gate voltage eV g /huJo = 8. Due to 
a finite £, the transition 0—^1 occurs at positive bias 
voltages. Moreover, the populations of the vibronic states 
are sketched in Figure fT5l for gate voltage V g = 0, which 
clearly shows that, for the considered parameters, only the 
vibronic ground state and first excited state are populated, 
whereas the populations of higher excited states are very 
small. This is in contrast to the non- stationary case where 
the excited states are highly populated (see Figure [5]). 
In a similar way, the populations of the vibronic states 
for gate voltage eV g /huJo = 8 are presented in figure [16] 
where higher excited states also get populated. Finally, 
in Figures [T7] and [18] we show the populations of the 0- 
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Fig. 13. (Color online). Population of (a) the 0-particle elec- 
tronic state, (b) the 1-particle electronic state. The value of 
gate voltage is V g = 0, and the frequency of the driving field 
is uj <C 1/Tmax- The other parameters are the same as used in 
Figure H 
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Fig. 15. (Color online). Populations as a function of normal- 
ized bias voltage for (a) vibronic ground state, (b) first excited 
state, (c) second excited state, and (d) fifth excited state. The 
frequency of the driving field is uj <C 1/ T max • The other param- 
eters are the same as used in Figure [2] 



Fig. 14. (Color online), (a) Population of the 0-particle elec- 
tronic state, (b) population of the 1-particle electronic state. 
The value of the gate voltage is eVg/hwo = 8, and the frequency 
of the driving field is uj <C l/r ma x- The other parameters are 
the same as used in Figure [2] 



and 1-particle vibronic states for gate voltages V g = 
and eVg/hujQ = 8, respectively, which basically provide 
the same information as mentioned before. 



7.1 I — V characteristics for the DC-case 

In the DC-case the analytical expression for the current 
remains the same as given by Eq. ([24]) taking into account 
a time independent bias. Let us first discuss the situa- 
tion when the 0- and 1-particle states are in resonance, 
e = eo — £ p = and V g = 0. In this particular case, an 
interesting behavior of the I — V characteristics with two 
opposite current peaks around zero bias can be observed 
(see Figure [19]) . In order to understand the mechanism of 
this process, we consider the source current which can be 
expressed in the form 

I B = T s F mm >{f + [e + fruo (m' - m) - eV b ] P m 

mm' 

-f-[e + huj (m f - m) - eV b ] Pf } • (42) 

At Vb = only ground to ground state transitions are open 
and P ° = = \. Hence, from Eq.(03J) one deduces that 
in this region the current is zero. At large positive bias, 
i.e., Vb — >• oo, the current is zero because the system is 
in a 1-particle stable state and no new transition channel 
is available. For finite bias, the behavior of the Franck- 
Condon factor F mm / is of importance. In particular, it 



0.95 



0.9 



0.85 



(a) ^ 






r 



q_ 0.05 




-60 -30 30 60 

eVb/Tiu>o 



0.015 



0.01 



0.005 



-60 -30 30 60 

eVb/fiooo 

x10" 3 





-60 -30 30 60 

eVb/fiooo 



-60 -30 30 60 

eVb/Tiooo 



Fig. 16. (Color online). Population as a function of normalized 
bias voltage for (a) vibronic ground state, (b) first excited state, 
(c) second excited state, and (d) fifth excited state. The value 
of the gate voltage is eV g /fwJo = 8, and the frequency of the 
driving field is uj <C 1 / r m ax • The other parameters are the same 
as used in Figure [2] 



suffices to investigate the classically allowed transitions as 
determined by the Franck-Condon parabola j32J:53J. The 

minimum of the parabola is for m = m' ~ (-|) , i.e., 

m,m' < (f) 2 transitions are exponentially suppressed. 
Moreover, F mrn ' attains the maximal values for F m7n > = 
F rn Q or F mrn > = Fom' and m or m! of the order of A 2 . 
Hence Fig. [T7] describes a threshold effect. The populations 
Pf 1 of the 1-particle states are mirror symmetric with re- 
spect to the bias inversion (not shown). Analogously, we 
can analyze in the same way as above the current peak 
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Fig. 17. (Color online). Populations P™ as a function of nor- 
malized bias voltage for the 0-particle (a) vibronic ground 
state, (b) first excited state, (c) second excited state, (d) third 
excited state, (e) fourth excited state, (f) fifth excited state, 
(g) sixth excited state, (h) seventh excited state, and (i) eighth 
excited state. The frequency of the driving field is uj <C l/r m ax- 
The other parameters are the same as used in Figure [2] 
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Fig. 18. (Color online). Population as a function of normal- 
ized bias voltage for the 0-particle (a) vibronic ground state, 
(b) first excited state, (c) fifth excited state and (d) 1-particle 
ground state, (e) first excited state, (f) fifth excited state. The 
value of the gate voltage is eV g /huJo = 8, and the frequency of 
the driving field is uj <C l/rmax. The other parameters are the 
same as used in Figure [2] 



in Figure [20] which occurs at eVh 
eV g /huJo = 8. 



e for gate voltage 



8 Conclusions 

In conclusion, we analyzed the quantum switching, bista- 
bility and memory effects in a single level system within 
the framework of the polaron model, where the electronic 
state is weakly coupled to metallic leads under AC-bias 
and strongly coupled to a vibrational mode. We showed 
that the bistability arises if the quantum switching be- 
tween neutral and charged states involved is suppressed, 
e.g., due to Franck-Condon blockade. In the case of an 
asymmetric junction, the neutral and charged states can 
be unstable at one polarity but stable at the other polarity 
of bias voltage. Under an appropriate choice of parame- 
ters, the stability regions of the two states overlap, which 




eV b /fi(j0 

Fig. 19. (Color online). Left current as a function of normal- 
ized bias. The frequency of the driving field is uj <C l/rmax- 
The other parameters are the same as used in Figure [2] 



x1CT 



1.5 



^, 1 



0.5 







J 





-60 -30 30 

eVb/ficjo 



60 



Fig. 20. (Color online). Current as a function of normalized 
bias. The value of gate voltage is, eV g /huJo = 8. The other 
parameters are the same as used in Figure [2] 



results in a bistable region in a certain interval of bias volt- 
age. Taking into account non-stationary effects, in particu- 
lar the interplay between time scales of variation of the ex- 
ternal perturbation and the switching time of the system, 
we demonstrated electrically controlled hysteretic behav- 
ior of the system. Furthermore, we showed that vibronic 
states and average vibron energies also show hysteretic 
behavior like the ones shown by the population-voltage 
and current-voltage curves. At the end, we also discussed 
the case of a DC-bias. In this case the population-voltage 
and current- volt age curves get single valued. Interestingly, 
one can observe current peaks in the I — V characteristics 
of the system when given vibronic channels contribute to 
transport. Moreover, we found that in the AC-case the 
vibronic excited states can be highly populated, while in 
the stationary case the population of the excited states is 
strongly decreased. 
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A Calculation of the correlator F a (t — t ; , /i ) 

Here we calculate the correlation function F a (t — t'j/xo) 
of lead a in the wide band limit. From Eq. (fT6|) we can 
write: 

F a (t ~ if, fJL ) =^2 fa (e K ~ floyW-V 

K, 

/co 
de£> a /(e-/i )e < *( t - t '> 
-CO 



= e *T?(*-*') 



dsD a f( 



if(i-t') 



(43) 



where D a is the constant density of states of lead a. To 
simplify the above equation, we use the following relation: 

/(g) = Hl+ e i"^ .Hence 



e 2 + e 2 



dee**'*"*') 



V 2 



(44) 



The first term leads to the result 



/CO 
dee »f(t-t') = 2irM(t-t'). 
-co 

The simplified form of the second part in ([44]) reads: 



(45) 



/^M -if (t-t') 



de tanh ( - — ) e l *> 



de tanh ( — — 
2 



x {cos +^sin ( 46 ) 

Due to symmetry the cosine component of the integral 
vanishes. One can further use the following relation: 



poo 

/ dx sin(ax) tanh 
Jo 



bx 



2 J 6tanh(^) 



, for a, b e R. 



Putting all together, the correlation function gets the final 
form (in the wide band limit) as 



F a (t - t'j/io) 



(49) 



h/3 



This function characterizes the correlation which exists on 
average between events where a lead electron is destroyed 
at time t' and another is created at time t. It thus provides 
very important information about the time scales which 
control the relaxation dynamics of the leads. 



B Evaluation of an integral 

In order to solve Eq. (|20|) and obtain the populations of 
the many-body states, one needs to evaluate the following 
integral: 

I a =J^dt"{F[t",^(t)]e-i AEt " 
+F*[t",i, a (t)]ei AEt "} 

= 2 J™ df'Re {F [t", fi a (t)] e -* AEt " } . (50) 

Substituting the correlation function F[t" ,/i a (i)] using 
Eq. (|4^|) in the above equation, we obtain 



J 



&t"-KhD a < COS 



-(n a (t)-AE) 



5(t") 



sin[^(^ a (t)-AE)] \ 
7i/3sinh(7r^) j 



(51) 



To simplify the above relation, one can use the following 
formula: 

sm(ax) 7rtanh(^) 
dx . ,\, = r-^-, a,b e M,and b > 0. 



j — c 



smh(bx) 



Using Eq. (|52)h we can write Eq. ([51]) as 



I a = 7rhD a < 1 + tanh 



13'- 



After some calculations, we obtain 

I a =2irHD a f[AE-ti a (t)], 
where jj, a {t) = /x + A/x a (t). 



(52) 
(53) 
(54) 



(47) 



Using (JUj), one can evaluate (|^6|) to be 



, ( P e 
detanh ( — — ) sin 



:(*"0 



2tt 



f3 sinh 



hp 



(48) 



C Evaluation of transition matrix elements 
of the electron operator 

To determine the transition rates, we need to calculate the 
matrix elements 

.1 I \|2 . 



(r\d\s) =e-2 |A| F(A,m,m'). 



(55) 
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where |r) and \s) represent the eigenstates given by Eq. (|9]). 
The function F(A, m, m') determines the coupling between 
states with a different vibronic number of excitations with 
effective coupling A and is expressed as [42 ,54J 

F(A, m, m ) = 

[e(m ; - m)A m - m + 9(m - m') (-A*) m ~ m '] 

(-iAi 2 r 



rn 



min • \ ^ 



+ m max - m min )! (m min - i)\ ' 



(56) 



where m min / max = min/max(ra, ra'). The coefficient F mm / 
inEqs. ([27j) and ([28]) is defined as F mm / = e~ x2 F 2 (\,m,m f ), 



D Expression for the FC factor F 2m f 

Using Appendix [C] the expression for F2m / follows to be 



A 



-A 4 + 2A 2 



\2(m'-2) 

+ ^r[ m ' 2 - m ' (1+2A2) + A4 ]^- (57) 
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